// # periodogram
// Оценка спектральной плотности мощности периодограммой
//
// В этом пакете мы объявляем функцию periodogram,
// которая принимает на вход сигнал x в виде среза []float64,
// окно window в виде среза []float64 и количество точек nfft для быстрого преобразования Фурье.
// Функция возвращает периодограмму в виде среза []float64.
//
// Мы начинаем с создания среза pxx для хранения периодограммы итерационно перебираем сигнал x в блоках размера nfft.
// Для каждого из этих блоков мы применяем окно window покомпонентно, умножая его на элементы блока x.
// Затем мы вычисляем быстрое преобразование Фурье (FFT) для окончательного преобразования блока x в частотный домен.
// Мы вычисляем квадрат модуля каждого элемента FFT и добавляем его к соответствующему элементу pxx.
//
// Также мы объявляем функцию fft, которая рекурсивно вычисляет быстрое преобразование Фурье для сигнала x.
// Эта функция работает за O(n log n) и использует формулу Баттерворта-Тьюки.
// Мы начинаем с проверки базового случая, когда n равно 1, и возвращаем значение x в качестве комплексного числа.
// В противном случае мы разделяем входной сигнал на две половины even и odd и рекурсивно вычисляем FFT для каждой половины.
// Затем мы объединяем два FFT с помощью преобразования Баттерворта-Тьюки и возвращаем результат.
package periodogram

import (
	"math"
	"math/cmplx"

	"github.com/mjibson/go-dsp/fft"
)

// Расчёт модифицированной периодограммы с окном
func Periodogram(x []float64, window []float64, nfft int) []float64 {
	pxx := make([]float64, nfft/2+1)
	for i := 0; i <= len(x)-nfft; i += nfft {
		block := x[i : i+nfft-1]

		//fmt.Println(i, len(x)-nfft, len(block), len(window))
		for j := range block {
			block[j] *= window[j]
		}
		FFT := fft.FFTReal(block)

		// fmt.Println(len(FFT), FFT)

		for j := range pxx {
			abs := cmplx.Abs(FFT[j])
			pxx[j] += math.Pow(abs, 2) / float64(nfft)
		}
	}
	return pxx
}
